/// Based on work by Thomas E. Kurek - C# WGS-84 implementation of zip code radius search
/// http://www.codeproject.com/KB/dotnet/Zip_code_radius_search.aspx

using System;

namespace ScrappyDB.Utilities
{
    public static class LatLonMath
    {
        //Equatorial radius of the earth from WGS 84 in meters, semi major axis = A
        private const int A = 6378137;

        //flattening = 1/298.257223563 = 0.0033528106647474805
        //first eccentricity squared = E = (2-flattening)*flattening
        private const double E = 0.0066943799901413165;

        //Miles to Meters conversion factor (take inverse for opposite)
        private const double MilesToMeters = 1609.347;

        //Degrees to Radians converstion factor (take inverse for opposite)
        private const double DegreesToRadians = Math.PI/180;

        public static ProximityBox CalculateProximity(double lat, double lon, double searchRadius, bool isMetric)
        {
            var result = new ProximityBox();

            //convert miles to meters
            if (isMetric == false)
            {
                searchRadius = searchRadius*MilesToMeters;
            }
            else
            {
                searchRadius = searchRadius*1000;
            }

            //lat naught and lon naught are the geodetic parameters in radians
            double lat0 = lat*DegreesToRadians;
            double lon0 = lon*DegreesToRadians;

            //Find reference ellipsoid radii
            double rm = CalculateMeridionalRadiusOfCurvature(lat0);
            double rpv = CalculateRadiusPrimeVertical(lat0);

            //Throw exception if search radius is greater than 1/4 of globe and thus breaks accuracy of model (mostly pertinent for russia, alaska, canada, peru, etc.)
            if (rpv*Math.Cos(lat0)*Math.PI/2 < searchRadius)
            {
                throw new ApplicationException(
                    "Search radius was too great to achieve an accurate result with this model.");
            }

            //Find boundaries for curvilinear plane that encloses search ellipse
            result.LatMax = (searchRadius/rm + lat0)/DegreesToRadians;
            result.LonMax = (searchRadius/(rpv*Math.Cos(lat0)) + lon0)/DegreesToRadians;
            result.LatMin = (lat0 - searchRadius/rm)/DegreesToRadians;
            result.LonMin = (lon0 - searchRadius/(rpv*Math.Cos(lat0)))/DegreesToRadians;

            return result;
        }

        /// <summary>
        /// Calculates the Radius of curvature in the prime vertical for the reference ellipsoid
        /// </summary>
        /// <remarks>
        /// This is the vector that defines the normal surface to any point on the ellipsoid.  It extends from
        /// from the polar axis to that point.  It is used for the longitude, in differentiation of east distances, dE
        /// </remarks>
        /// <param name="lat0">Geodetic latitude in radians</param>
        /// <returns>Length of radius of curvature in the prime vertical</returns>
        private static double CalculateRadiusPrimeVertical(double lat0)
        {
            double rpv = A/Math.Sqrt(1 - E*Math.Pow(Math.Sin(lat0), 2));
            return rpv;
        }

        /// <summary>
        /// Calculates the meridional radius of curvature for the reference ellipsoid
        /// </summary>
        /// <remarks>
        /// This is the radius of a circle that fits the earth curvature in the meridian at the latitude chosen.
        /// It is used for latitude, in differentiation of north distances, dN
        /// </remarks>
        /// <param name="lat0">Geodetic latitude in radians</param>
        /// <returns>Length of meridional radius of curvature</returns>
        private static double CalculateMeridionalRadiusOfCurvature(double lat0)
        {
            double rm = A*(1 - E)/Math.Pow(1 - E*(Math.Pow(Math.Sin(lat0), 2)), 1.5);
            return rm;
        }

        /// <summary>
        /// Calculates the true birds-eye view length of the arc between two positions on the globe using parameters from WGS 84,
        /// used in aviation and GPS.
        /// </summary>
        /// <remarks>
        /// An accurate BIRDS EYE numerical approximation with error approaching less than 10 feet at a 50 miles search radius
        /// and only 60 ft at 400 mile search radius.  Error is on the order of (searchRadius/equatorialRadius)^2.  Only accurate for distances
        /// less than 1/4 of the globe (~10,000 km at the equator, and approaching 0 km at the poles).
        /// Geoid height above sea level is assumed to be zero, which is the only deviation from GPS, and another source of error.
        /// </remarks>
        /// <param name="Rm">Meridional Radius of Curvature at the centroid latitude</param>
        /// <param name="Rpv">Radius of Curvature in the Prime Vertical at the centroid latitude</param>
        /// <param name="lat0">Centroid latitude</param>
        /// <param name="lon0">Centroid longitude</param>
        /// <param name="lat">Destination latitude</param>
        /// <param name="lon">Destination longitude</param>
        /// <param name="isMetric">true for kilometers, false for miles</param>
        /// <returns>Distance in miles or kilometers from the arc between the two points on the globe</returns>
        public static double CalculateDistance(double lat0, double lon0, double lat, double lon, bool isMetric)
        {
            lat0 = lat0*DegreesToRadians;
            lon0 = lon0*DegreesToRadians;
            lat = lat*DegreesToRadians;
            lon = lon*DegreesToRadians;

            double rm = CalculateMeridionalRadiusOfCurvature(lat0);
            double rpv = CalculateRadiusPrimeVertical(lat0);

            double distance =
                Math.Sqrt(Math.Pow(rm, 2)*Math.Pow(lat - lat0, 2) +
                          Math.Pow(rpv, 2)*Math.Pow(Math.Cos(lat0), 2)*Math.Pow(lon - lon0, 2));

            if (isMetric == false)
            {
                distance = distance/MilesToMeters;
            }
            else
            {
                distance = distance/1000;
            }

            return distance;
        }
    }

    public class ProximityBox
    {
        public double LatMax { get; set; }
        public double LonMax { get; set; }
        public double LatMin { get; set; }
        public double LonMin { get; set; }
    }
}